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Introduction / 3 * 

In a recent series of studies [1-3J it was shown that the 
method of expansion in terms of spherical harmonics is a 
powerful tool which can be applied to solving the transfer 
equation in scattering media with scattering characteristics 
varying over the optical depth. Jeans [^1] proposed the idea 
of this method; it was further developed during its appliea- 
tion to the problem of the transfer of neutrons (corresponding 
bibliographical references are given in [1]), Canosa and 
Penafiel [1] applied the method to the transfer of radiation 
in anisotrotropically scattering homogeneous media and they 
pointed out an approach by which cumulative rounding errors 
can be eliminated during its application to media with great 
optical depth. Dave and Canosa [2] demonstrated the applica- 
bility of the method to inhomogeneous media. All these studies 
considered only the azumith--independent radiation intensity 
component (zero-th azimuth harmonic). Subsequently Dave [33 
elaborated a generalization of the method which can he used 
to find higher azimuth harmonics. 

This study has two objectives. First, all studies per- 
taining to the discussed method with which the author is 
familiar are. based on using so-called, .’.’.vacuum. 5 ', .bou nda ry. 
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conditions on the upper and lower boundary of the scattering 
medium. When these conditions are applied to the atmosphere 
of a planet , they imply that the albedo of the lower boundary 
surface is assumed to be zero. Such a constraint is not jus- 
tified. Therefore relations which allow us to apply the given 
method to more realistic lower boundaries are derived in this 
paper . 

Second, all studies based on an application of the given /4_ 
method known to the author pertain to transfer of solar radia- 
tion, i.e. the assumption is made that the source function is 
determined hy radiation from a source of parallel rays located 
outside the scattering medium. In this paper the case of 
thermal radiation transfer is also considered, i.e, it is 
assumed that the source function varies with depth and that it 
depends on the temperature of the atmosphere and the emissive 
power of the source, 

I. Expansion in Spherical Harmonics 

We present here without derivation the basic relations 
used: in this method applied to the case of solar radiation 
transfer with ’’vacuum" boundary conditions. We follow studies 
[1, 2] and retain most of the notation used in them. Differences 
pertain only to a number of minor, misprints and also to some 
changes made for more convenient electronic computer programming. 

Let us consider a scattering and absorbing medium having 
optical depth whose scattering and absorption characteristics 
depend on the depth. These characteristics can be described 
by specifying the coefficients to^, in the expansion 

of the scattering indieatrix in terms of Legendre polynomials 

P ( cos 0) = Z. P, C cos 
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normalized in such a way that the quantity ^ represents the ' 
single scattering albedo. The selected number L is odd [l3‘. 
The original transfer equation together with the ’’vacuum" 
boundary conditions is expressed in the form: 


< 2 


dl ° ^ 

or. 


I(T. ; U/'f ) = 

-2PF j \ uift u',f) I(r ; u'jf') du'&f' 

^ „ * 

* ? P ®2Cp(~ ~ ) PC*^ 5 VL » f } ^ } 

* 

I( 0? uy ^ ) = 0 when ~ u > 0 
I<r g? uj<f ) s 0 when 


u < O 


/5 

( 1 ) 


where I is the intensity of the scattered radiation; 

t is the optical thickness (depth) measured from the 
upper boundary of the atmosphere; 

Uj u’ is the cosine of the zenith direction angle of 
scattered radiation measured from the positive direction 
of the t axis (nadir); 

P is the scattering indicatr-ix; 

irP is the illumination intensity of an area beyond the 
atmosphere which is perpendicular to the direction of 
the Sun’s rays; 

P 0 is the cosine of the Sun’s zenith angle measured 
from the zenith; 

$0 is the azimuth of the Sun. 

The azimuth independent .intensity component 

. "i * y 

I (“? u) J I (v: Uifj df (2) 

o 
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is sought in the form of an expansion in Legendre polynomials 
(henceforth the superscript n 0 ,f will be omitted) 


where 


I(x- 5 u) = J - ) P-^(u) j 

ho 


m 

f 1 ( t ) = ^ I c $ u) P- OO 3» t 




and the problem is reduced to fin ing the vector 


r t f-, O 


f(X)r 




i 

tr(~ ) 


f iC-) 


( 3 ) 

( 4 ) 


Equation (1) is reduced to a system of differential equations 
with boundary conditions, which in matrix notation, has the c '/§_ 
form: * * ^ 

A * c.(*& )f( - ) = s(^). 


<l~o 

f o (0) = -Gf e (0), .. 
f o (^).x.Gf e (xp, 


( 5 ) 


where A and C(t) are matrices having rank (L+l): 


± - 


A 

0 ^ 

^2> 0 A 


0 


L-1 


3 # 

0 


.0 


( 6 ) 
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whose elements are 

, _ 1 1 . p 1 .. 

~ -h = ~ 5 ~2T-r\" ? 

1 ' » 


C7) 


C8) 

(9) 


respectively 3 and s(x) is an L+l-dimensional vector with compo- 
nents 


S e (~ 


) = 


• F - 


Cl0) 


fg and f^ are (L+l)/2— dimensional vectors, whose components 
are respectively the odd and even components of the vector f: 


✓ > 

- 


£, 1 

- 


?3 

« 

o *P * * — 

» r e “ 

?2 

a 

% 



< -J 

J 

k. ^ 


G is the matrix of boundary conditions found on the basis of /7 
the recurrence formulas 


G 1,0 

G 1 ? 21 



(41*1X21-3) - - 




1-2 


5 J>5 


*■ © o o o o 




. (ID 


0 21 z 




r 2m-1 ? 21 


J- 


= 2$ 3g oci »«© (Iri-i)/2 „ 

Elements G^ with even i and odd 4 are not defined, and the 
matrix is considered to have rank (L+l)/2. 
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The next step is the introduction of a grid .of n integra- 
tion points T^ s with t^O and T n “ T g» an< ^ reduction of (5) to 
a system of matrix equations having the form 


+ S i 0 i-c.1 f i-r1 “ w i ’ 


2. = 


_ A 


© & b 


n-1. 


(f 1 ) o = - 


"Ke s ^n-'o v “*n'e > 


( 12 ) 


where 



S i 9 iv1 = A * 




r— — — °C (— — — ) t 


w. = 
1 


~' z i 


£ s ( J . 


(13) 


Next , row and column operators are introduced, which interchange 
rows or columns of the matrix according to‘ eveness , i.e. in 
such?: a way, that the first (L+l)/2 rows (columns) of the matrix 
are even (0, 2, L-l) rows (columns) and the following 

(L+D/2 rows (columns) are odd (1, 3,..., L) rows (columns)* 

The system of matrix equations (12) is transformed to the form 



•°i 9 i s i * 

/8 

( 14 ) 


i- “ ^ y 2 $ oooo 5 X ^*1 ^ 



ev= - GSi 5 s£ = GSa, • 

( 15 ) 

where 

D 1 ^ » Column (S -^) 9 

( 16 ) 



= How (£ ± V 
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The superscripts ”t" and t'b 1 ' denote respectively the top or 
bottom half of a vector or matrix. The initial system of 
equations for obtaining a solution is system (l 4 3 15 ), 

Canosa and Penafiel [ 1 H presented methods for solving 
equations [14, 15 ] both for the case of moderate optical 
depth and for media whose optical depth is so great that it 
requires certain stabilizing transformations preventing an 
increase in rounding errors. 


2.1, Medium with Moderate Optical Depth 

A sequence is constructed from n-1 matrices having rank 
(L+l): 




1 = In 2 a 000 ■ 


original; page 1$ 

OS SOOR QUALITY 


( 17 ) 


A sequence is constructed from n rectangular (L+l) x (L+l)/2 
matrices 






i = 2j> 3s 


>n 


( 18 ) 


where. I. is the identity matrix. 

A sequence is constructed from n (L+l) -dimensional vectors: 






0 

L 


"1 = +■ 

^ 3 $ 000 n » 


( 19 ) 


The (L+l )/2-dimensional vector^. is determined from a solution 
of the system of CL+l )/2 equations having the following form 
in matrix notation: 
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The sequence consisting of the n (L+l) -dimensional vectors 


g- - ir 1 * a, 

— ^ i 


T — * 

UK' "" t J 


‘X > 

2 < i c- O O. o J 


Y9 


( 21 ) 


is a solution of system (14, 15). 


2,2. Medium with Great Optical Depth 

The necessary number of so-called, overdetermining points 
ti (k=l , 2.. , N c ) with i =1, ij| =n, is selected from the ' 

set of points x£. The points are selected on the basis of the 
requirement that the rounding errors cannot increase to an 
appreciable magnitude on an interval between neighboring over- 
determining points. Empirical selection rules for the position 
of the points are discussed in [2]. 


A sequence is constructed from N transformation matrices 

“ V 

having rank (L+l)/2: 





■li * 

k = 2g J? I 



( 22 ) 


A sequence is constructed from n transformed (L+l)/2— dimensional 
vectors 


4- o 

8 


= a i, *1 ^1 2 2°2 * 


K 


o p « **- T/^2° ° ^ 

k ss 3s A-p ooo » • / • 


(23) 


A sequence is 
matrices : 


constructed from n rectangular (L+l) x (L+l)/2 


u i = 3?,i « o 

i<:L k . 

e o c^L 


fc=2 s 3, 



= W-2 


T 

« c o o o— ° 


(18-1) 
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A sequence is constructed from n (L+l ) -dimensional vectors: 

f 0 ^ 


, a i^ * *** * ^T 2CSO T fe t k ), 




'• i"* 2 V 3 


h 


$ o * # 


SL 
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The (L+l )/2-dimensional vector 1 m is determined from a solution 

c 

of the system of (L+l)/2 equations having the following form 
in matrix notation: 


ci& 


«2»* 



$ 20 - 1 ) 


The sequence consisting of N — 1 (Lt-1 )/2=dimensional vectors 1, : 

c K 

* . ‘ # " • # * - 

^k “• ®k+1 ^k+1 


k.= F c -1» 3ff 0 -2 9 o 


(20-2) 


is determined 

The sequence of n (L+l)— dimensional vectors g^ representing 
a solution of system (l4, 15) is determined: 


Si - ^i^k * u i i \ 1 

^ ^k+1 » J f 

k B 1, 2 9 o.c (21-1) 


• s a a ’ a n m 

3. Boundary Condition Matrix for Lambertian' Lower Boundary 

Formulas for elements of the matrix of boundary conditions 
(11) were presented above for the case when no scattered radia- 
tion impinges on the scattering medium from outside. When the 
medium is bounded (for example from below) by a boundary with 
nonzero albedo 3 such radiation is obviously present. Let us 
assume that the optical depth of the scattering atmosphere is 
sufficiently great and the illumination intensity of the lower 
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boundary due to direct solar radiation- is negligibly small. 

Then the boundary condition for the Lambertian boundary with 

albedo A can be written in the form: 

4' 2# /ii 

, x ( 1 W f .) = | \ \ I(-| g u-f fj 

0 0 

This relation can be transformed to the condition for the azimuths- 

independent component expressed in the form: 

• * 

s 2& J •. u p ) -d*dtt‘° 

' ' ' ' 0 - ** ' - ’ - 
taking into account (3), we have 

• • \ - l ' ' ' * 

u'au' = $ %<•%) p x ( u5Kdu'= 

'0 o 

. - L 1 ^ 

, =Z [ i *,( «') «') te'] W -Z°u£[<V 


(24) 


Formulas for the C ml are given, for instance, in [2], For odd 
1, the following holds: 


°na = - 1 


0,5 .when !a = 1 

0 .when 1 


(25A) 


For even 1, the quantities C jnl and the elements of matrix G 
defined by formulas (11) are related by the simple relation 

Gal = 2C ml (25B) 


Following the scheme proposed by Marschak (see for instance 
[2]), writing down expansion (3) for the left member of (24) 
and taking into account subsequent mathematical operations we 
obtain : 


L 


z. 


21+1- 




P,( u ) = 24. z: 

. ■ .. i-0 . , ■ 


( 26 ) 
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$*t/Gi2v A . u 

Since (26), can oiitly be satisfied for a finite number of points 
in [-1, 0 J , following Marschak, we multiply both sides of (26) 
b y (m=l, 3,..., L) , integrate over the interval [-1, 03, ,4 ' 

and require that these equalities be satisfied for all m. We 
have ovfl- • 

£ l(z . g) Pl (u) P B (u) du = 2 A \ 4a. 

ift.fter obvious transformations, writing 


: D ml. T 




we obtain: 




m = i 9 3s 


(27) 


Performing obvious operations it can be shown that 


A for 


even 


C ml A for odd 


Taking into account (25A), we rewrite (24) in the form 


1 .* 

■ .>2 A a- 


i' ■ Wi = -" AC m o 4 £ i '■ 

?=0 2 ... * „ , T 


•ms 1 . 3 


Let us first consider the case ra=l. Substituting C^q= 0.25, 
after obvious transformations and taking into account (25B), 
we obtain : 


( 28 ) 


-P - 

1 " 1?K 




G 11 f l 


(29) 
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Next, let us consider the case of arbitrary odd -m. Substituting 
(23) in (28) and taking into account (26) werobtain: 



(30) 


It can be easily shown that (29) is a special case of (30) for 
m=l . Relation (30) holds for m=l, 3, . .., L, and it can be 
rewritten in matrix formr 


• = %f e (-=g) C3D 

where the elements of matrix G^ have the form: 

^A^ml = G ml ~ T53T %bO G 11 ‘ 

We obtained the sought form of the boundary condition for 
the Lambertian boundary surface* The matrix must be sub- 
stituted in equations (20) and (21) instead of the matrix G. 

It is easily seen that 


G g+G as A-»-0 

Another obvious property follows: when A=l, having performed 

elementary operations, we have: 

CG A )qi*=° C33) 

This corresponds to the fact that at a boundary with A=l, we A3 
must have zero radiation flow. Indeed, it can be easily seen 
that the magnitude of flow at the lower boundary is airf^Cr ). 
Combining (31) and (33)* we obtain 

fi (T g ) - 0 

i,e. the flow is indeed zero. 

^ Thermal Radiation Transfer 

The thermal radiation transfer equation can be written in 
the following form (see for instance [53): 

114 



0E FOOBl 


u* 


dlC-c 5 u) 


* i \ P( T I u * u ') I(T * 

^ -1 


+ I(-C 5 U) « 

u') du' ‘-4- [l- <%«>] B C^) 


( 34 ) 


where B(T) is the radiation intensity of an ideal black body. 
As before, the upper boundary condition is a vacuum condition, 
and the lower boundary condition can be written in the form 

"T 

ICv = 2A I I<V U*) u'du' + (1-A)B(!T S ) 

o 

where 1-A, T"- is the emissive power and temperature of the 
' s 

lower boundary respectively. 


Next, we follow= the general idea of the spherical harmonics 
method. Applying the summation theorem for Legendre polynomials 
(see for instance [63) and changing to the variable t in the 
function B(T), we rewrite (34) in the form: 

u * lC*c ? u) = 

= ^ P 1 (u) ^ 1^(0 5 u')dtt'+ 

Multiplying the above equation by P k (u) (k=0, 1, ..., L) and 
integrating from —1 to 1, we obtain after obvious transforma- 
tions 


A - 1 

} f§ p k (u) tt du + [ 1 “ !£&£•] 5 u5 p k Cu)du 


-1 

k = 0, 1, »o. L 


where 


/14 


b k C^) 


, [l - a 0 (t) ] • B(x) 

* 

0 


for k = 0 

for k = i j 2 9 o.oI i- 
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Applying the recurrence relation 


(2k+1) u P k (u) = (k4-1) P k ^(u) + kP k _^(u) 


(see for instance [?]) and replacing the subscript k by 1 3 we 
obtain the system of equations 


The validity of the system written in this form can be verified 
by keeping in mind that P^=0 f° r a -^ an< ^ taking f^^O- 
Here 3^ and o^(t) are determined on the basis of formulas 

(8, 9), and for s^Ct) werhave: 


+ <31 = 8 1 C c ') 


1 = 0 , 


(36) 


*!<■*> = 2151 * 


2[l - ai t Ct)].B(T) 


' l for 


System of (36) can be written in matrix notation in a form which 
is identical with the form of the differential equation in 
system (5). Thus, as was to be expected, the differenc from 
the case of solar radiation transfer is the source function 
vector s(t) which is now determined from formula (37) instead 
of (10). 

Let us now consider the lower boundary condition (35). 

The right member involves a constant which is independent of 
I(t ; u) within the context of the given formulation of the 

O r 

problem. Denoting for conciseness this constant °y b, and pro- 
ceeding analologously as in section 3, we write: 
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L we have ; 


For m«l, 3, l 


f— D ml f l " 4AD mO ' C 1i f l + 

6 L-1 


L-1 


If ~ y C = -4AC ~(4 c^f,) 

2 H t%,~ ml 1 102 1 f -'A- 11 1 

<®AA- 


2C md b 


V&0 y 2 ... 


=RE- G mO b 


From the above we obtain the sought expression for the boundary 
condition in matrix form: 

•W = ®aVV - tIt 

Here G Q is the i-th column of the matrix G. 

Clearly, as A-^0 and B(T_)->0, (38) becomes the "vacuum” 

boundary condition for t in system (5). 

S 

5. Conclusions ' 

It was shown above that the spherical harmonics method 
can be generalized to the case of a lower boundary with albedo 
different from zero and that it can also be applied to the 
case of thermal radiation transfer. In both cases the effect 
of the lower boundary can be taken into account by modifying 
the lower boundary condition matrix. In addition, in the case 
of thermal radiation transfer, it is necessary to make an 
appropriate replacement of the source function vector. The 
above-mentioned changes are minor and they can be carried out 
as additional modifications of existing computer programs using 
the spherical harmonics method. 

The author is deeply grateful to V.I. Moroz, doctor of 
physical and mathematical sciences, and V.G. Zolotukhin, 
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